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Recently it was shown [A. Gordon and F. X. Kartner, Phys. Rev. Lett. 95, 223901 (2005)] that 
the strong field approximation (SFA) for high-order harmonic generation (HHG) is significantly 
improved when the SFA wave function is used with the acceleration rather than the length form of 
the dipole operator. In this work it is shown that using the acceleration form upgrades the SFA 
from zeroth-order to first-order accuracy in the binding potential. The first-order correct three-step 
model (l st -order TSM) obtained thereby is systematically compared to its standard zeroth-order 
counterpart (0 th -order TSM) and it is found that they differ significantly even for energetic electrons. 
For molecules (in the single-electron approximation), the th -order and the l st -order TSMs in general 
disagree about the connection between the orbital symmetry and the positions of the minima in the 
HHG spectrum. At last, we briefly comment on gauge and translation invariance issues of the SFA. 

PACS numbers: 32.80.Rm, 42.65.Ky 



X 



I. INTRODUCTION 

The strong field approximation (SFA) 0- 13 is a 
key technique in the study of interactions of matter with 
intense laser fields. In particular, the SFA is used to 
describe high-order harmonic generation (HHG) 0, 0, 0, 
%, The three-step model (TSM) [HEQ, which is based 
on the SFA, has proven very successful describing much of 
the experimental behavior. The best known examples are 
the cutoff formula [RJp, Q and time-frequency structure 
of the HHG signal ||, which led to the prediction of 
attosecond pulses @- Reviews can be found in Refs. 0, 

It has been long known however that the TSM is incor- 
rect by 1-2 orders of magnitude predicting the spectral 
intensity of HHG in atomic hydrogen, as found from com- 
parisons with numerically-exact results 0, ^2] • For the 
H^" ion the TSM gives a spectrum which is many orders 
of magnitude away from numerically-exact results |13| . 
and the shape of the spectrum can be heavily distorted 

mm. 

HHG experiments are now coming to the point where 
more accurate theory is needed. The race towards har- 
HHG for a coherent short wavelength source 
can benefit from quantitative theoretical esti- 
mates of the HHG efficiency, which to date are very scarce 
[l8l | . The recent orbital imaging experiment 0] uses the 
HHG spectrum to infer the structure of a molecular or- 
bital, which also requires quantitatively reliable theory, 
capable of giving a precise description of the shape of the 
HHG spectrum. 

In a recent theoretical work [l^ | we proposed a modi- 
fied version of the TSM, where the SFA wavefunction is 
used with the dipole operator in the acceleration rather 
than the length form. Comparison of the TSM with a 
numerical solution of the time dependent Schrddinger 



equation (TDSE) then demonstrated excellent quanti- 
tative agreement for atomic hydrogen, and significantly 
improved agreement for the H^" ion. The argument why 
using the acceleration form improves the TSM so much 
was that the TSM then becomes correct to first order 
in the binding potential (we shall henceforth refer to the 
TSM obtained this way as l st -order TSM), whereas the 
standard TSM HH0 (O th -order TSM) is correct to ze- 
roth order in the binding potential. 

This work is a followup on Ref. an d ^ nas three 
goals. First, the first-order accuracy of the l st -order TSM 
is established. Second, a detailed and general comparison 
between the l st -order TSM and the O th -order TSM is 
given. It is shown that the two always disagree with 
respect to the connection between the orbital symmetry 
and the positions of the minima in the HHG spectrum, 
as demonstrated for H J [TliL IbH ] . Third, the opportunity 
of re-deriving the SFA is used to give a perspective on 
some delicate issues in the derivation, such as gauge and 
translation invariance, which have been under debate ppl 

mm mm. 

This paper is organized as follows. In Sections ILT1 and 
ITTT1 the derivation of the SFA and th -order TSM is re- 
viewed, in an attempt to illuminate some subtleties in 
the derivation and in order to set the stage for deriving 
the l st -order TSM in Sec. |TV| Section |V| compares the 
O th -order and l st -order TSMs, and SecEJis dedicated to 
the discussion of gauge and translation invariance issues 
in the SFA. Section IVlII gives a brief summary. 

Atomic units are adopted throughout the paper. 



II. STRONG FIELD APPROXIMATION 

This work is restricted to the single-electron approxi- 
mation, where an atom or a molecule are modeled by an 
electron in an effective (local) potential V(r): 
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K(r) 



(1) 
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I p is the binding energy of the ground state, and is added 
to Eq. for convenience reasons, such that the ground 
state of Ho has zero energy. The atom or molecule is 
placed in a linearly polarized electric field E(t). The 
x axis is chosen along the direction of polarization, the 
wavelength is assumed sufficiently long such that the 
dipole approximation holds, and the length gauge is cho- 
sen, to give the Hamiltonian 



H(t) =H Q - E(t)a 



(2) 



The SFA is usually presented [see, e. g. Ref. 0] and 
references therein] as a perturbative expansion in V(r), 
where the unperturbed Hamiltonian is the Volkov Hamil- 
tonian 



Hv(t) = --V 2 -E{t)x + I p . 



(3) 



The evolution operator Uv(t,t') of Hy(t), defined by 

id t U v (t,t') = H v (t)U v (t,t'); U v (t',t') = 1, (4) 

is known exactly |6J. Using the Lippmann-Schwinger 
equation one can approximate U(t,t'), the evolution op- 
erator associated with H (t) (defined through Eq. Q with 
all V subscripts omitted) to arbitrary order in V(r) pfij . 
The zeroth and first order would be 



U {t,t') = U v {t,t') 



(5) 



U x {t,t') = -i / dt"U v (t,t")VU v (t" ,t'). (6) 



With the operator U(t,t') at hand, any problem can be 
solved. 

If at time t = the electron is in the ground state |0) of 
H a , the wavefunction \ip(t)) at any time would be given 
by 



!</>(*)> = [/(i,0)|0>. 



(7) 



One therefore could suggest approximating U(t, t') in 
Eq. (J7J by the perturbative expansion outlined in Eqs. JS] 
|HJ). However this is not the way the SFA is usually per- 
formed. The route is rather by making the ansatz j(| 



M=a(f)|0> + |p(t)>. 



(8) 



a(t) has the initial conditions a(0) = 1 and is determined 
later. In other words, the way we split into a(t)\0) 
and \<p(t)) is not yet defined. 

Let now \<p(t)) be the exact solution of 

#(*)> = H(t)\<p(t)) - a(t)E(t)x\0), (9) 

which in terms of U reads 
t 

\<p(t)) =-i I U{t,t')E(t')xa(t')\0)dt'. (10) 



Let us assume for the moment that a(t) is specified. The 
SFA approximates U in Eq. IjlOjl by a perturbation series 
in V upon XJy. 

One can now ask why the SFA approximates U in 
Eq. (|10|l rather than in Eq. JjJ, if both lead to an exact 
wavefunction when U is exact? In other words, why one 
needs the ansatz (|SJ| and the rather non-straightforward 
Eq. iPJ? To our understanding, the answer to these ques- 
tions is that the perturbed Uv(t,t') is a very bad ap- 
proximation for describing the evolution of (quasi) bound 
states. 

Without the laser field, perturbation theory for the 
evolution operator (the Born series) diverges at bound 
states of the full Hamiltonian 26] . The presence of a laser 
field may change the situation, since formally there are 
no bound states at all. We are not in a position to make a 
rigorous mathematical statement about the convergence 
of the SFA, but since the ground state remains quasi- 
bound, one may expect convergence difficulties. Obvi- 
ously, at least for the the zeroth order theory, Uv(t,t') 
itself gives a very bad approximation for the evolution of 
a quasi-bound state. In fact, in the Keldysh ionization 
rate yj is found rather indirectly, by calculating the rate 
at which the norm of the continuum increases, and the 
rate at which the ground-state decays is inferred solely 
by conservation of the total norm. 

The last term in Eq. © is therefore added because we 
know in advance that we are going to use an approximate 
propagator, which will give a very bad description for the 
evolution of the ground state. Yet, Eq. provides a 
clear route for systematically improving the SFA: a(t) is 
first assumed to be known, and Eq. @ is solved to in 
principle an arbitrary order in V(r). Then a(t) is found 
by demanding the conservation of norm, or by means of 
other approximations, as the perturbed Uv(t,t') alone is 
not well suited for this purpose. Note that in the limit 
where U is exact, one easily finds from Eq. © and Eq. (JSJ) 
that a(t) = 1. 



III. ZEROTH ORDER TSM 

In this section we re-derive the standard th -order 
TSM. The main reason is that our slightly different 
derivation lays the technical foundations for establishing 
the first order accuracy of the l st -order TSM in Sec. IIVI 

The th -order TSM is obtained by simply replacing 
U(t,t') by U v (t,t') in Eq. ©. Since U v {t,t') is known 
exactly, the solution can be written in a closed form [(|: 

(p-A(t)l^(t)) = 
t 

= -i [ dt'a(t')E(t')(p - A(t')|x|0)e^ lS ( p ' t '*'). (11) 



A(t) = —x f E(t')dt' is a vector potential that describes 



3 



the electric field E(t), and 



z 

S(p,t,1?) = ± J(p-A(t")) 2 dt"+I p (t- 



t<). (12) 



At this point we slightly deviate from the standard 
derivation of the TSM 0, Q: We perform the saddle- 
point integration in Eq. (|llfl right now, before proceed- 
ing. This calculation was pioneered by Keldysh |lj and 
is further discussed in a vast number of SFA studies [see 
e. g. Ref. and references therein]. The discussion 
in this work is limited to the tunneling regime, which is 
defined by the requirements [27^ : 



E « (2J P ) 



W <C In 



3/2 



(13a) 
(13b) 
(13c) 



In Eq. I|13a|) E is the amplitude of the driving field and 
lo is its frequency. These parameters are precisely de- 
fined for a sinusoidal driving field, and more loosely for a 
general field, such as a sinusoidal field with an envelope. 
In the latter case, u> can be thought of as the parameter 
characterizing the timescale over which E(t) varies. 7 
is the well known Keldysh parameter p|. Note that the 
requirement (|13c(l follows from Ijl3a(l and l|13b(l . 

Under the conditions (|13|) , the saddle-point integration 
of Eq. Qll[) is carried out to give 



( P -A(tMt)) = ^53°(*»(P*)) 



X g \B(t n (p X ))\y 



J p e 



w((E(t n (p x )) 
\E(t n (p x ))\ 

-iS(p,t,t n (p x )) ^ 



w(E) is the static Stark ionization rate associated with 
the ground state, p\ = p 2 + pi, and the function t n (p x ) 
is the set of positive real solutions of the equation 

p x =A x (t n (p x )). (15) 

The number of solutions depends on E(t). a(t) can now 
be found by requiring conservation of the norm in Eq. © 
[neglecting (0|y(i))] to give the well-known expression 



|a(i)| 2 = e -/o 



(16) 



Using the ansatz (JSJ with the wavefunction (|14f) . one 
can now compute the expectation value of the dipole mo- 
ment: 

(x) = \a(t)\ 2 {0\x\0) + { V (t)\x\^(t))+^(t)+C (t),(17) 
where 

t o (t)=a*(t)(0\x\<p(t)). (18) 

The origin of coordinates can always be chosen such that 
the first term in Eq. (|17fl vanishes. The second term 



has no high harmonics 28], since matrix elements of x 
between two different Volkov states vanish. The high 
harmonics come from the cross term £0 (t) . Using Eq. i|14|) 
we find: 



&(t)=a*(t) 




d 3 p(0\x\p - A(t))a(t n ( Px ))x 



w((E(t n (p x )) — jg(p t itn(px)) , . 
\E(t n ( Px ))\ ^ 1 ' 

The p integration is now carried out in the stationary 
phase approximation as in Ref. Q . The stationary phase 
is attained at p± = and at all p x values such that t n (p x ) 
satisfy [see Appendix lAl for a more careful discussion] 



(A(t n ) - A(t"))dt" 



(20) 



For any given t, i n (t) is defined as the set of solutions 
to Eq. I|20(l . i. e. birth times of trajectories that end up 
at the origin at time t. Then in the stationary phase 
approximation Eq. I|19|) gives 



iS n (t) 



{t-t n (t)f/-- 



;a*(t)x 



x a(t n (t))(0\x\A(t n (t)) Mt)) W ^l^\ (21) 

where S n (t) = S(A(t n (t)), t, i n (t)). 

Eq. (|21l) is identical to the original TSM expression de- 
rived in Ref. [7| . The only difference is that Eq. (|21|l takes 
into account the depletion of the ground state, and that 
all numerical prefactors (2 3 / 2 7r) were calculated. The re- 
derivation of the th -order TSM is now concluded. 



IV. FIRST ORDER TSM 

A straight forward way to upgrade the TSM from ze- 
roth to first order accuracy in V would be improving the 
TSM wavefunction (|llf> by adding the first-order correc- 
tion to the evolution operator, Eq. (JSJ. This method has 
been employed in Ref. '29], for improving the SFA the- 
oretical description of above threshold ionization. For 
HHG we propose another method, which is significantly 
simpler, and does not require correcting the wavefunc- 
tion. 

If \if>(t)) is the exact solution of the time-dependent 
Schrodinger equation with the Hamiltonian J2J , then the 
Ehrenfest theorem holds: 



dt 2 

d(m\px\m) _ 

dt 

-(m\d x V(r)\m)+E(t) 



(22a) 

(22b) 
(22c) 
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Therefore computing the time-dependent dipole expec- 
tation value in all three forms of Eq. I122II . length (|22a[) . 
velocity |j22bj, or acceleration (j22c|l . is equivalent if ac- 
companied by the appropriate differentiation or integra- 
tion in time. 

However with the approximate wavefunction given by 
Eqs. © and (|14|) . the results will be in general different 
for each form. One therefore has to make a choice which 
form to use for this particular wavefunction. We argue 
that the acceleration form gives in general the best re- 
sults, since even if \f(t)) is only correct to zeroth order 
in V(r), the expectation value is automatically correct, 
at least formally, to first order in V(r). One therefore 
needs to evaluate 



dt 2 



£i(«) + £(*)+£c(t) 



where 



&(t) = -a*(t)(0\d x V(r)\v(t)) 

Ut) = -{ V {t)\d x v(v)\ v {t)) 



(23) 



(24) 
(25) 



The last term of Eq. (|22c|l has been dropped, since it is 
the driving field itself and thus contains no harmonics. 
It is easy to show, exploiting the commutator [p, Hq], 
that (0|9 K y(r)|0) = [irrespectively of the position of 
origin] , which is why the corresponding term in Eq. I|23|) 
is missing. 

Note that the mere appearance of V in Eq. (|23|) is 
not sufficient to warrant first order accuracy. One can 
see that, for example, from the fact that zeroth order 
SFA expressions in V can be transformed such that they 
appear linear or quadratic in V [23> . The (formal) 
first-order accuracy of Eq. (|23|l is established by the fact 
that a first-order correction in V to the wavefunction in 
Eq. will result in only a second order correction in V 
in Eq. (|23|) . This is why we dedicated Sec.[n]to carefully 
defining the procedure by which the SFA is corrected 
order by order in V. 

The cross term £i has the same structure as £o 
[Eq. (|18[l ]. Going through the same steps that have led 
to Eq. I|21(l . one arrives at the same expression £i, with 
the only difference that the matrix element (0|x|k) is re- 
placed by -(0\8 x V(r)\k): 



&(*) 



v 1 4" (t-t n (t))v* 



x a* (t)a(t n (t)) 



w((E(i n (t)) 



\E(in(t))\ 
x (0\d x V(r)\A(t n (t))-A(t)). 



(26) 



Eq. I|26|) is the improved version of the TSM presented 
in Ref. The change in the expression compared to 

Eq. H21JI is very small and easy to implement, and yet 
results in a very large difference, especially in the case of 
molecules. 

Eq. 126(1 is an expression for £i, but Eq. Il2.'ill contains 
also £ c . It turns out that the contribution of £ c to HHG is 



smaller by at least 0(w 3//2 ) than that of £i. £ c is therefore 
negligible for w<l, which coincides with Eq. I|13cf) when 
I p is of 0(1). The latter holds for all neutral (or not 
highly charged) atoms and molecules. The evaluation of 
£ c is rather lengthy, especially for potentials with a long- 
ranged Coulomb tail, and is given in Appendix iBl 

V. th VS. l st -ORDER TSM - COMPARISON 

The O th -order TSM suggests Eq. i|18|) as an approxima- 
tion to the dipole moment, whereas the l st -order TSM 
suggests Eq. (|24H . In order to compare them conve- 
niently, we now differentiate Eq. I|18|) twice in time and 
compare £o with £i. This is done in detail in Appendix 
lO The result is that under the condition (|13afl . in or- 
der to obtain £o(t), one nas to replace the l st -order TSM 
recombination amplitude 



<= c w (k) = (0\d x V(r)\k), 



in Eq. (gljl by 



,old 



(k) = (/ p + I|k| 2 ) 2 



(27) 



(28) 



Comparing the l st -order TSM and the O th -order TSM is 
thus reduced to comparing af™ and a* respectively. 

The two expressions look different, and indeed they 
are. A detailed comparison requires the knowledge of 
V. However in order to gain some general insight, in 
what follows we study the k — > oo asymptotic behavior 
of a"™ and a*. Some pretty general statements can be 
made about k — > oo, which turn out to provide important 
insights also for k-s of 0(1). 



A. High-momentum asymptotic evaluation 

a°l^(k) can be written as 



<c(k) 



1 



X old (r) e - lkr dr 



(2tt) 3 / 2 

(and similarly for a" c c w (k) and x ncw (r)), where 



(29) 



X ° ld (r) = (/ p -iV 2 ) 2 ^ (r) (30) 



X new (r) = Mr)d x V(r), 



(31) 



with ?Ao(r) = (r|0). Using ff |0) = 0, Eq. gJJ can be 
transformed into 



X old (r) - 2Mr)dxV(r) + 2V(r)d x Mr) 

+ xV 2 (r)M*) + zW(r) • V^o(r)+ 
+ ix^o(r)V 2 ^(r) 



(32) 



The k — * oo asymptotic behavior of a°l^.(k) and a"™(k) 
is dictated by the most singular part X° ld ( r ) an d X ncw ( r ) 
respectively, since they are connected by a Fourier trans- 
form m. 
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We now assume that V(r) is an effective potential that 
represents a molecule (or, as a special case, an atom) with 
nuclei at positions r, with charges Zj, j — 1...N. Com- 
mon choices for an effective potential to model such a 
system have a Zj/\r — Tj\ singularity near the j-th nu- 
cleus, which is partially screened off away from the nu- 
cleus. Then the last term in Eq. 1|32|1 can be written as 

ixV>o(r)V 2 y(r) = -27ra^o(r) £ Z 3 S{v - Tj ) (33) 

3 

Now one has to distinguish two cases. The first case is 
a molecule. Whenever there is more than one nucleus, 
Eq. (|33|l does not vanish [apart from possibly a very few 
orientations of the molecule] . Then the term in Eq. (|33|l is 
the most singular one among all terms on the right hand 
side of Eq. {jJSJ, and thus governs a°^(k) at k — > oo. In 
this case 

X old (r) ~ -2irY,XjMrj)Zj5(r - Tj ) (34) 

3 

where the notation '~' means 'equal in its most singular 
part'. 

The second case is an atom, which has only one nu- 
cleus, whose position can be always chosen at r = [see 
Sec. lVll for a discussion about translation invariance]. In 
this case Eq. (I33|) obviously vanishes. Then the most sin- 
gular term on the right hand side of Eq. \.\'2\) is the first 
one, and we find that 

X old (r) ~ 2 X ncw (r) (35) 




X 



FIG. 1: comparison of arc™ [thick red lines] and a"^ [thin 
blue lines] for hydrogen. The continuous lines correspond to 
the two expressions in Eq. (OHJ. The dashed lines are the 
asymptotic expressions [Eq. 1361 and Eq. 1371 1. The inset 
shows the ratio between the two expressions in Eq. Q380 . 



C. Molecules 



B. Atoms 

Eq. (|35|l implies that for atoms 

< d (k) - 2< c w (k), (36) 

where '— >' denotes 'approaches asymptotically at k — > 
oo'. Thus a"°"(k) and a°' d (k) do not agree even asymp- 
totically for k — > oo, the limit were plane waves approx- 
imate the exact continuum states increasingly well. For 
atoms, the O th -order TSM predicts a 4 times larger HHG 
yield than the (l st -order TSM). For a Z/r singularity in 
the potential, the asymptotic behavior of a^™ can be 
easily worked out to give 

a^(k)^iZM0)^l^ (37) 

When -00(0) vanishes Eq. (|37|) is modified. 

For finite k the disagreement between a" c c c w (k) and 
a°,l d (k) can be even greater. For the hydrogen atom one 
finds |32l ] 

new/ 7 \ k x tan k x 



The differences between the -order TSM and the 1 st - 
order TSM for molecules is particularly striking. The two 
models give opposite relations between orbital symmetry 
and the positions of the minima in the HHG spectrum, 
a subject of many recent theoretical and experimental 
studies [H HI HE IH HI- A set of minima 

is predicted by the th -order TSM with an odd orbital, 
corresponds to an even orbital in the l st -order TSM, and 
vice versa. 

From Eq. (|34|l we find that for molecules 

a°H(k) --4=E xjM^Zse-^, (39) 
V j 

whereas from Eq. I|27|l and Eq. I|37|) we find 

<7(k) ^iyf^E^(ri)^e- lk -i, (40) 

3 

Observing Eq. l|3^|l and Eq. l|4Tj|l . it becomes clear why 
the th -order TSM and l st -order TSM disagree about the 
connection between the symmetry ?/>o and the zeros (or 
minima) of |a rec |. Apart from the additional k x /k 2 en- 
velope in Eq. I|4l)[). which does not affect the position of 
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the zeros, Eq. (|40|l is the Fourier transform of xipo(r), 
whereas Eq. i|39|) is the Fourier transform of ipo(r) (both 
sampled at the singular points). xtpo{ r ) and ^o( r ) ob- 
viously have the opposite symmetry with respect to a 
x — > —x reflection. 

In order to illustrate this we consider the H^ ion, where 
the nuclei are positioned along the x axis, at x = ±R/2. 
In this case Eq. gives 

a^(k)^^sin(±Mi), (41) 

where C = i()o(—^Rx) — ipo(^Rx). In contrast, Eq. l|4TJ|) 
gives 

a^(k)^i^^2Ccos^k x R), (42) 

A zero of a rcc at a given k x corresponds to a minimum 
in the HHG spectral intensity at the frequency I p + 
The O th -order TSM therefore predicts minima at energies 
given by I p + (2irn/R) 2 (n is an integer), whereas the 
l st -order TSM predicts minima at energies given by I p + 
(2it(n + i)/i?) 2 . The latter condition well agrees with 
numerically-exact results 0, |3j| . 

It is easy to see from Eq. l{3*9")) and Eq. (|4Tj)l that if ip (r) 
were an odd rather than even wavefunction with respect 
to x — > —x, the cosine [sine] in Eq. I|41|l [Eq. (|42[l ] is re- 
placed by a sine [cosine 1. Therefore if the O th -order TSM 
is used for reconstructing ipo(r) from the HHG spectrum, 
it flips the symmetry of tp (r) from even to odd and vice 
versa. This statement is, of course, based only on the 
k x — ► oo asymptotic behavior of the a rcc -s. Yet it is 
interesting to note that for Hj, Eq. (|4*T)l and Eq. J32J| 
approximate Eq. (|27|l and Eq. I|28(l reasonably well fol- 
low momenta, as Fig. [21 shows. 
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FIG. 2: \a°f c \ 2 (thin blue line) and |a° e c w | 2 (thick red line) 
as function of v x for (R=2). The dashed lines are the 
corresponding asymptotic (v x 3> 1) approximations [Eq. 14111 
and Eq. gg)]. 

Note that the k x /k 2 envelope, which is present in 
Eq. (|40|l but absent in Eq. (|39f) . results in an orders- 
of-magnitude difference in the magnitudes of the a rcc -s, 
in addition to the distorted shape. 



VI. GAUGE AND TRANSLATION 
INVARIANCE ISSUES 

The dependence of the SFA on gauge has been the 
subject of lively discussions H |U HH H i3| . Here 
we do not attempt a comprehensive discussion of this 
subject. However since some steps in Sec. IV! seem to rely 
on placing the origin of the laser potential and the atom 
at the same point, we briefly address the related issue of 
translation invariance. 

Consider thus the Hamiltonian 

H(t) = -\v 2 + V(r) - E(t){x - x ) + I p (43) 

(xq is a constant), which is related to Eq. (J2J) by a gauge 
transformation. The exact solution of the TDSE with 
Eq. obviously gives the same time-dependent expec- 
tation values as with Eq. J5J). However the SFA can give 
different results for Eq. l|43fl and Eq. @. In particular, 
the SFA with Eq. 143fl can give rise to the generation of 
even harmonics [21|, which is an artifact. 

There is an obvious way by which this problem can be 
cured, namely, by replacing the ansatz (jSJ by 

\xp) = a(t)e~ lx ° f E ^ dt \0) + \(p(t)). (44) 

Substituted Eq. JUJ) into the TDSE with the Hamilto- 
nian (|43|) , it is easy to convince oneself that the results 
no longer depend on Xq. This is because wavefunction 
(JSJ) undergoes exactly the gauge transformation that con- 
nects Eq. (|43|) to Eq. @, guaranteeing gauge invariance. 
The SFA depends on the choice of gauge simply because 
it crucially depends on the choice of the initial ansatz. 

One can now ask how to obtain an ansatz for a given 
gauge a-priori, without e. g. gauge-transforming Eq. JSJ, 
and how to obtain the ansatz (JSJ) itself a-priori. Our an- 
swer is that the ansatz should provide the best approxi- 
mation for the evolution of the quasi-bound ground state, 
since the (perturbed) Volkov evolution operator cannot 
do it [see Sec. HT]. 

In the low frequency regime [Eq. (|13c|) ]. a good ap- 
proximation for the evolution of the ground state can 
be found using the adiabatic theorem |4j : It is approx- 
imately given by < E W) dt \ <j>(E)}, where the \<j>{E)) 
and s(E) are the field-dependent ground state wavefunc- 
tion and energy, defined through the eigenvalue equation 

(_ i V 2 + V(r) - E(x - x ) + I p ) \4(E)) = e{E)\<j>{E)), 

(45) 

The adiabatic theorem holds even though \<fi(E)) be- 
comes a resonance K we neglect the Stark shift 
of the energy, we obtain e(E) = Exq. Going further and 
approximating \<f>(E)) by |^(0)) = |0), we obtain exactly 
the ansatz (|44|l . a-priori, without referring to Eq. JSJ). In 
fact, Eq. (JSJ itself is obtained for xq = 0. 

The discussion can be extended to at least one common 
gauge, namely the velocity gauge. Then the adiabatic ar- 
gument would again give a modified ansatz, similar to the 
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one used in Ref. |4^. It is likely that the differences be- 
tween the length and velocity gauges reported in Ref. j2^| 
would then disappear. 



VII. DISCUSSION 

The large discrepancy between the O th -order TSM 
HHG spectra and numerically-exact calculations has 
been attributed long ago to the SFA wavefunction 
being inaccurate. It is often argued that this SFA wave- 
function is especially inaccurate for molecules |4fll |. This 
work strongly supports these ideas: The length and ac- 
celeration forms of the dipole operator used on the same 
SFA wavefunction are shown to give dramatically differ- 
ent results, especially for molecules. For an exact wave- 
function the results would be identical, and thus the large 
discrepancy is evidence for the inaccuracy of the wave- 
function. 

Since the wavefunction is so far from being exact, it 
is especially important to use it properly. In this work 
we have argued that together with the acceleration form 
of the dipole, the zeroth-order SFA wavefunction leads 
to an approximation (the l st -order TSM) with first-order 
overall accuracy in the binding potential. Since the wave- 
function itself is so inaccurate, it can lead to very large 
errors if used otherwise than in the specially-suited way 
introduced in this work. The accuracy of the l st -order 
TSM should be tested in terms of the expectation values 
it generates and not in terms of the wavefunction itself. 

We have shown that the continuum-continuum term in 
the l st -order TSM is negligible compared to the bound- 
continuum term. Therefore, the widely-used O th -order 
TSM is upgraded to first-order accuracy by simply re- 
placing x by d x V(r) in the recombination amplitude, 
which is simple to implement. The l st -order TSM shows 
excellent agreement with numerically-exact results for 
atoms 0, EH and good agreement for . 

The SFA decomposes the wavefunction into a contin- 
uum, and a quasi-bound ground state. This is justified, 
since the strong laser field smears the excited states of the 
unperturbed system and turns them into a continuum . 
The dynamics of the continuum is approximated by the 
Volkov propagator, which can be corrected order by order 
in the binding potential. The quasi-bound state however 
should be treated separately. 

It should be noted that there could be more than one 
quasi-bound state. For example, the initial state of the 
electron can be an excited state of the laser-free Hamil- 
tonian a case that was not treated in this work. This 
is commonly the case in single-electron models of multi- 
electron atoms. Then the states lying energetically be- 
low the initial state are also quasi-bound, and should 
also be treated separately. One simple way to do that is 
to project these states out of the SFA wavefunction, as 
briefly discussed at the end of Ref. 0] . 
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APPENDIX A: COMPLEMENTARY REMARKS 
ON THE DERIVATION OF THE th -ORDER TSM 

Eq. I|2(J|) is gives the stationary phase condition only 
to leading order in 7. In fact, when S(p,t,t n (px)) is dif- 
ferentiated with respect p x taking Eq. I|15|) into account, 
the result is 

t 

^ry + J(A(t n )-A(t"))dt" = 0. (Al) 

tn 

The second term in Eq. (|Alj) is of the order of Eo/u) 2 , 
while the first one is of the order of I p /Eq. Therefore 
the first term in Eq. (|A1I) is a second-order correction in 
7 to Eq. H2Ufl . It is interesting to note however that the 
subleading term term is meaningful, since the electron is 
released by tunneling roughly at the turning point {x = 
Ip/E) [27| rather than at the origin. 

Note that the denominator (t - t n (i)) 3 / 2 in Eq. J5lJ 
cannot result in a divergence. If the first term in Eq. 12U|) 
is retained, the reason is obvious - the duration of a tra- 
jectory is never zero, since it begins and ends at different 
locations. If the first term in Eq. (|20|l is neglected, as 
often happens, one can easily show that trajectories for 
which t approaches t n (t) can only occur at points where 
E(t) — 0, and then w(E) vanishes exponentially. Per- 
forming the t' integration first thus spares the need for 
the regularizing parameter e used in, e. g. Ref. |f|. 

Note also that the denominator (t — t n (t)) 3 ^ 2 is cor- 
rect only to zeroth order in 7. The leading correction is 
obtained if the e^^^^P^V 2 ^ factor in Eq. O is 
included in the Gaussian in integration and if Eq. (|15f) is 
taken into account when the second derivative of S with 
respect to p x is computed. 



APPENDIX B: EVALUATION OF THE 
CONTINUUM-CONTINUUM TERM 

For completeness, in this Appendix the continuum- 
continuum term £ c (i) [Eq. 1|25[) ] is evaluated. To this 
aim we define 

Q(k) = J d 3 r[d x V(r)]e- ik *, (Bl) 

the Fourier transform of the force field derived from the 
potential V. It turns out that a long-range Coulomb be- 
havior of V requires more careful evaluation of £ c than a 
short-ranged V. Since once \(fi(t)) is given £ c (t) is linear 
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in V , one can write V as a sum of the pure Coulomb po- 
tential plus a short-ranged potential, evaluate each term 
separately, and sum up the results at the end. In what 
follows we therefore treat the two cases separately. 



2. Coulomb potential — different trajectories 

Now we consider the case V(r) = — 1/r, which leads to 



1. Short-range potentials 

In order to keep the expressions from becoming too 
cluttered, we introduce the notation 



a( f ) = __ a( i)__ 



(B2) 



Substituting Eq. lfTI|) in Eq. (|2l>jl we obtain 

L(t) = J dpdp'a*(t n ,(p' x ))a(t n (p x ))x 



l-B(t„(p? c ))l % /2T 



p e 



p x 



Q(p' - p) (B3) 



We now carry out the integration in Eq. I|B3|) in the 
stationary-phase approximation. To this aim we assume 
that apart from the exponentials in Eq. (|B3(1 . the rest of 
the integrand is slowly varying in p. Here it is where we 
use the short-range property of V, which assures that Q 
is indeed slowly varying in p. 

The stationary phase condition in the six-dimensional 
p-p' momentum space is identical to the one of the inte- 
gral {T5)| : p± = p' ± = 0, whereas p x and p' x are obtained 
by finding all solutions i n (t) that satisfy Eq. 12U|) and 
using Eq. (fT3f) to find the corresponding p x -s. The result 
of the integration is 



a* (t n ,)a(t n ) 



(t-f„)3/2(t_ t - n/ )3/2 



xe 



iS(t,t n ,)-iS{t,t n ) 



Q(A(t„0-A(t n ))), (B4) 



where the dependence of t n on t has been suppressed. 

Eq. 1B4|) confirms our claim that £ c (i) is small com- 
pared to £,.(*). [Eq. ffift] has only a (27r/(t - 
tn)) 3 ^ 2 prefactor (which is 0(w 3 / 2 )), whereas for £ c (t) 
it is 0(lu 3 ). Moreover, £i(t) is proportional to y / w(E), 

whereas t; c (t) is linear in w(E). HHG experiments typi- 
cally operate under the condition of small ionization per 
cycle, which means w(E) <C ui. Therefore although £ c (i) 
is formally of the same order in V(r) as £ c (i) is 

more than 0(ui 2 ) smaller than and is thus negligi- 

ble under the conditions 

It should be noted that since is proportional to 
a* (t), at very high fields, where the ground state is almost 
completely ionized in one cycle, £i(i) becomes exponen- 
tially small in the field amplitude 0]. I* 1 this case £ c (i) 
may become significant. However this regime is of little 
interest form the point of view of HHG, since HHG ba- 
sically disappears under these operating conditions |44| . 



Q(k) = i 



7T k 2 



(B5) 



Due to the singularity at k = 0, the stationary phase 
approximation should be now used more carefully when 
integrating Eq. (|B3|) . To this aim we now consider an 
integral of the form 



Ip-pT' 

(B6) 

/ represents the third line of Eq. (|B3|) . and the rest of 
the integrand is slowly varying and can be added later. 
Eq. (|B6|) has two oscillating phase factors, centered at 
po and Pq, and the Coulomb potential, r and r' are 
real and positive, and by comparison with Eq. I|B3(I one 
can see that they represent the traveling times of the two 
trajectories. 

Using the convolution and Plancharel's theorems, 
Eq. 1B6(1 is transformed to 



where 



(2ir) 3/2 - 

1 = ^im I ^ s ^ 



o xexp(i^ iApox) 

a r . 



(B7) 



(B8) 



where s = 1 /t' — 1/r and Apo = po x — p' 0x . We assumed 
that the po — p[) is parallel to the x axis, since this is the 
case of interest in Eq. (|B3(I . 

The integral in Eq. i|B8|) can be carried out analytically 
in spherical coordinates and expressed in a closed form 
using the error function: 



47ri 




Ap Q . (B9) 



Eq. l|B9Jl and Eq. I|B7J| give an exact expression for / in 
Eq. I|B6|) . Let us now perform the integration in Eq. JB6I) 
using the stationary-phase approximation instead. The 
result is 



-^SPA = 



i2 7 / 2 7T 5 / 2 

Apr 3 /V 3 / 2 ' 



(BIO) 



Since each exponential picks up only an environment of 
radius r -1 / 2 (or r' -1 / 2 ) around its center po (or p^), one 
expects that when 1/r, 1/r' <C Ap 2 ,, the integration does 
not reach the singularity and Eq. (|B10|) holds. Eq. I|B7|) 
and Eq. i|B9(l verify this expectation. Figure visualizes 
this statement, and shows that I(Ap 0l s) basically gives a 
smoothed version of the 1/Apo singularity near Apo = 0. 
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The discussion in Sec. IB ll therefore holds as it is as long 
as the condition 1/r, 1/V <C Ap^ is met. This condition 
is violated when Ap approaches zero, and this case will 
be our concern in what follows. By observing Eq. I|B4|) 
one can see that the latter always happens when n = n' , 
and can also accidentally happen if n 7^ n! . We begin 
with the second case. 



stationary phase approximation and can be avoided if one 
uses the fact that the first two lines in Eq. (|B3|I vanish 
exponentially at p,p' — > 00. Doing so is however techni- 
cally cumbersome, and we adopt a simpler approach. 

We start with having another look at Eq. (|B4|I . For a 
given pair n and n' , the main HHG frequency to be gen- 
erated is given by the derivative of the phase in Eq. (|B4|I : 




FIG. 3: |/(Ap ,s)| [Eq. (Eg)] as function of Ap - The thick 
line is 47r/Apo, the limit of |J(Apo, s)| at s = 0. 



By Taylor-expanding the modulus of Eq. (|B9J| keeping 
the two leading orders in Apo, one can find the maximal 
value of |J(Apoj s)| with respect to Apo, and see that it 
is proportional to s -1 / 2 . Using this result for an upper 
bound on /, one obtains 



\I\ < 2 7 ' 2 



700 



r 5/2 



80919 • 



(Bll) 



Note that we are considering two different trajectories, 
which means that r ^ t' . Moreover, it is easy to show 
that two different trajectories that return at the same 
time must have their birth times separated by more than 
a quarter of a driving cycle [u>\t — r'| > n/2}. It follows 
therefore that Eq. l|BlT|l is 0(w 5/2 ). 

For two different trajectories, n ^ n' , the Coulomb 
singularity therefore enhances the integral Eq. (|B3jl by 
at most a factor of uj^ 1 / 2 . This is seen by comparing 
Eq. (|B7|) and Eq. (fBTTfl . In Sec.|B7]we have shown that 
£ c (t) is smaller than by a factor of 0(uj 2 ). Now we 
arrive at the conclusion that in the long range case, for 
two different trajectories, £ c (^) can be sometimes smaller 
than £i(i) by factor of 0{lo 3 / 2 ) only. Yet, £ c (i) remains 
negligible for u> <C 1. 



3. Coulomb potential — same trajectory 

If n = n', which means r — r' and Ap Q = 0, the inte- 
gral l|B6|) diverges. This divergence is an artifact of the 



dt 



(S(t,t n ,)- 
= \{A{t) 



S(t,t n )) : 

-A(t n ,)) 2 



- 2 (A(t) 



A{t n )) 2 (B12) 



Eq. (|B12|) has a simple intuitive meaning. gives the 
beat frequency between the ground state (at frequency 
Ip) and a continuum electron with a frequency corre- 
sponding to the kinetic energy upon return. In contrast, 
£ c (i) is obtained from two different electron trajectories, 
which return to the parent ion at the same time with 
two different kinetic energies. The emitted radiation is 
at the beat frequency corresponding to the difference in 
the kinetic energies upon return. 

Obviously, if n = n' , that is, if the two trajectories orig- 
inated at the same birth time, they are identical. There- 
fore the beat frequency is zero and no high harmonics 
are generated. In other words, for n = n' Eq. (|B3|) varies 
slowly in time. There is yet a subtle point to be checked: 
loosely speaking, if Eq. I|B3|I is very large, then even if it 
varies very slowly in time, its high-frequency tail may be 
comparable to £1 (t) . 

In order to show that this is not the case, we introduce 
an infrared cutoff to the Coulomb potential, replacing 
it by the Yukawa potential V(r) — — e~ r l r ° jr. After 
changing the variables of integration to p + = i(p + p') 
and p_ = p — p', Eq. (|B3|) has the form 



d 3 p+d 3 p-e" p+p 



. P-x9 '(P+iP- 
P- + ro 2 



(B13) 



where g represents the slowly-varying part of the inte- 
grand. 

The integral (|B13I) is governed by the vicinity of p+ = 
p_ = 0, and we therefore Taylor-expand g around the 
origin. The zcroth order term, where g(p+,p_) is re- 
placed by <7(0,0), vanishes since the integrand is odd in 
(p+jP-)- Since the first two lines of the integrand in 
Eq. i|B4|) consist of two identical functions of p and p', 
the derivative of g with respect to p_ vanishes at the ori- 
gin. The leading non-vanishing term in Eq. (|B13|) comes 
from the derivative of g with respect to p+ x , giving: 



iy -g + 



P-xP+x 



P- 



9+ 
3r 2 



d 3 p, 



3 -rp + /r 
P+ 



^g+rl 
3r 4 ' 



(B14) 



where g+ = 9 p+ ,<?(0,0) = -3(f„(i)) / 'E(t n (t)). 

The range of the Yukawa potential, ro, should be set to 
a value that is large enough such that the wavefunction 
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is contained within. Classically, the electron travels to 
distances of the order of E /uj 2 , whereas r is of the order 
of l/w. Therefore the modulus of expression in Eq. (|B14|I 
is of the order of a(i n (t))E . 

We shall not delve into the analysis of a(i n (t))Eo any 
further, being content with the fact that it is much 
smaller than one, and is a smooth function of time that 
varies over a timescale of u> . These two properties as- 
sure that frequencies of O(l) have amplitudes which are 
exponentially small in u , and thus can be neglected. 
Our statement, that the continuum-continuum term £ c (t) 
is negligible compared to the recombination term £1 (t) as 
far as HHG is concerned, is now established. 

APPENDIX C: DERIVATION OF EQ. lE8l) 

The zeroth order SFA wavefunction is defined through 
the equation 

i\<p(t)) = Hv(t)\tp(t)) ~ a(t)E(t)x\0). (CI) 

We now differentiate £o twice in time, using Eq. IjClfl . 
We obtain: 

&(*) = iEa*(t)(0\x 2 \^) - a*{t){Q\xH v {t) 2 \v), (C2) 

where terms including derivatives of a(t) have been 
dropped, as they are negligible. Terms that do not con- 



tain \ip) have been dropped as well, as they cannot con- 
tribute to HHG. 

If we substitute Eq. I|14|) into Eq. (|C2|) . we obtain an 
expression for £o(t) which is identical to Eq. H26[) . with 
the only difference that the matrix element (0|9 l: V r (r)|p) 
is replaced by 

(0\xH 2 (t)\p) - iE(t)(0\x 2 \p) = 
= + (0\x\p)-2E(t)(^+I p y0\x 2 \p) 
- E 2 (t)(0\x 3 \p)+E(t)i(Q\xp x \p)-iE(t)(0\x 2 \p) (C3) 

Of the terms on the right hand side of Eq. I|C3(1 . we ar- 
gue that the first one is the most dominant in the regime 
defined by Eq. (|13a[) . In order to show this we note that 
(0|p) has a pole at p 2 = —2I p . This follows simply from 
the asymptotic behavior of (0|r) at r — > oo, which is 
an exponential falloff at the rate of \J21 P . Therefore 
for p 2 > 2I p , (0|p) typically falls off like p 2 + 2I p to 
some negative power. This fact can be used to estimate 
the modulus of the second term on the right hand side 
of Eq. (EU as \E(t)y/p^+2I p '{0\x\p)\ up to a factor of 
order one. Using <C (2/ p ) 3 / 2 from Eq. (|13a|l one ar- 
rives at the conclusion that the second term is negligible 
compared to the first one. Using similar argumentation 
and using Eq. I|13afl one shows that the first term on the 
right hand side of Eq. (|C3|) is indeed the dominant one. 
Eq. is thus established. 
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